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Abstract. The Dirac equation is an important model in relativistic quantum mechanics. In the 
semi-classical regime e<^l, even a spatially spectrally accurate time splitting method [6] requires 
the mesh size to be O(e), which makes the direct simulation extremely expensive. In this paper, we 
present the Gaussian beam method for the Dirac equation. With the help of an eigenvalue decom- 
position, the Gaussian beams can be independently evolved along each eigenspace and summed to 
construct an approximate solution of the Dirac equation. Moreover, the proposed Eulerian Gaussian 
beam keeps the advantages of constructing the Hessian matrices by simply using level set functions' 
derivatives. Finally, several numerical examples show the efficiency and accuracy of the method. 
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1. Introduction We are interested in developing the Gaussian beam method 
for the Dirac equation in the semi-classical regime 

ied t ^ e = -tea ■ V* e - a ■ + /3* £ + VV E , (1.1) 

subject to the Cauchy initial data 

& e (0,x) = Ui(x) iSl & )/e , xgM 3 . (1.2) 

Here * fe (t,a;) = (*f,*|,*|,*|) !r eC 4 is the spinor field, normalized s.t., 

/ |* e (t,x)| a dx = l, 

0<£<1 denotes the semi-classical parameter, V(t,x) GM. is the external electric 
potential and A(t,x)GM. 3 represents the external magnetic potential, i.e., A = 
(Ai,A2, As). Without loss of generality, we only consider static external field in this 
paper. The Dirac matrices /3,a= (ai,a2,cxs) are complex-valued Hermitian matrices, 
which are given by 
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Gaussian beam for the Dirac equation 



with I 2 the 2x2 identity matrix and a k the 2x2 Pauli matrices, i.e., 

H")' M-)< •*-(;-.)■ 

The physical observables can be defined in terms of *& E (t,x): 

Particle-density p e = |* £ | 2 , (1.3) 
Current-density f k = (* e ,a' c * £ ) C 4 . (1.4) 

The Dirac equation [2J [35] is a relativistic wave equation which plays a funda- 
mental role in relativistic quantum mechanics. It provides a natural description of 
the particles with spin 1/2 , i.e. electrons, neutrinos, muons, protons, neutrons, etc. 
The Dirac equation also predicts some peculiar effects, such as Klein's paradox [13] 
and "Zitterbewegung" an unexpected quivering motion of a free relativistic quantum 
particle [29] . Recently, graphene [U [24] and topological insulators [5j [38] are studied 
widely in connection to the Dirac equation in the semiclassical regime }22| . 

In the semi-classical regime eCl, the solution to the Dirac equation is highly 
oscillatory. Thus, for any domain-based discretization method, the number of mesh 
points in each spatial direction should be at least 0(e~ 1 ) [14]. If the potential is 
sufficiently smooth, and the initial data of the Dirac equation is compactly supported, 
the time-splitting spectral method 6 offers the best numerical resolution. The spatial 
meshing strategy is almost of optimal order 0(e~ 1 ) and the time step can be O(l). 

One alternative efficient numerical approach for solving the Dirac equation is the 
WKB method [301 [HD EH] • For a first order approximation, this method tries to seek 
an asymptotic solution: 

^ e (t,x) = u(t,x)e iS( - t ' se ^ e + 0(e), SeR, (1.5) 

where the amplitude u and the phase S are smooth functions independent of e. Substi- 
tuting (|1.5[) into the Dirac equation, one derives the eikonal equation and the transport 
equation. Since the eikonal equation is of the Hamilton- Jacobi type, the solution be- 
comes singular after caustic formulation. Beyond caustics, the correct semi-classical 
solution of the Dirac equation contains several phases. In the last decades, many 
approaches have been proposed to capture this multi-phased solutions, see reviews 

A serious drawback of the WKB method is that the solution ceases to be valid 
at caustics where the rays intersect and the amplitudes blow up. The Gaussian beam 
method, which was first proposed by Heller in quantum chemistry 4 and indepen- 
dently developed by Popov in Geophysics [26], is an efficient approach that allows 
accurate computation of the amplitude and phase information near caustics. The 
main difference between the WKB method and the Gaussian beam method is that 
the Gaussian beam method allows the phase function to be complex off the center 
of the beam and the imaginary part of the phase function is positive definite, which 
makes the solution decay exponentially away from the center. The validity of the 
Gaussian beam method at caustics was analyzed by Ralston in [28]. The Gaussian 
beam and related methods have become very popular in high frequency waves prob- 
lems [111 [I9l [201 [231 [271 Ell Ell [391 [40] in recent years. Most of the methods were in the 
Lagrangian framework. More recently, Eulerian Gaussian beam methods have also re- 
ceived a special attention for its advantages of uniform accuracy [TU1 [TT1 [T2l [TBI [TBI [T7] . 
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A major simplification of the Eulerian Gaussian beam method is that the Hessian ma- 
trices can be constructed by taking derivatives of the level set functions [TU]. This 
greatly reduces the computational cost. 

To our knowledge, no Gaussian beam methods have previously been developed 
for the Dirac equation. It is the goal of this paper to develop such a method by 
extending the previous method of [lOj for the Schrodinger equation to the Dirac 
equation (|1.1[) - (|1.2[) . With the help of the eigenvalue decomposition, the Gaussian 
beams evolve independently of each other. Moreover, the energy transition is for- 
bidden since the Dirac matrix /? results an 0(1) band gap. Being different from the 
Gaussian beam methods for the Schrodinger equation, the higher order Taylor expan- 
sion and asymptotic expansion must be considered for the amplitude. After making 
use of the solvability condition and matching the different expansions, one gets the 
transport equation for the lower order term of the amplitude. When the evolution 
is done, the solution can be simply constructed by the summation of all Gaussian 
beams. The solution will be shown to have a good accuracy even around caustics, 
with a coarse mesh size of 0(y / e) and large time step of 0(y / e). A remarkable aspect 
of the Eulerian Gaussian beam method is that it still possesses the advantage of the 
previous method |10j . which is an important benefit for the ZD simulation. 

The paper is organized as follows. After reviewing the semi-classical limit of the 
Dirac equation in Section [2j we formulate the Gaussian beam method for (|l.ip - (ll.2|) 
in Section [3] In Section |4l our method is shown to be accurate and efficient by several 
numerical examples. Finally, we conclude the paper in Section [SJ 

2. The Dirac equation and the semi-classical limit Denote the Dirac 
operator by 

V(x,Z) = a-(Z-A(x))+P + V(x), 

then we have 

V (x,-ieV) * e = a (-ieV - A{x)) * e + /3* e + V{x)^ e . 
Therefore the semi-classically scaled Dirac equation (jl.l[) can be written as 

ied t ^ £ = V (x, -ieV) * e . 

Let 

\(x } $)=\f\a-A(x)\ 2 +i, 

then 

h ± (x,£) = ±\(x,£) + V(x) 

are two different eigenvalues, each with multiplicity two, of the Dirac operator T>(x,£). 
The corresponding projectors n ± (a;,^) are given by 

II^*,*) = \ (l± ± {V(x,Z) - V(x)h)j . 

Plugging the following WKB-ansatz into (|l.ip , 

oo 

* £ (t , x) = e iS ^/ e £ J ''Uj (* , x) , 



4 Gaussian beam for the Dirac equation 

where Uj(t,x) £ C°°(1R 4 ,C 4 ), matching the O(l) and 0(e) asymptotic coefficients, one 
has 

(d t S+V{x,VS))u = 0, (2.1) 
i(dt + a-V)u o -(dtS + V{x,VS))u 1 =0. (2.2) 

In order to get a nontrivial solution uo(t,x) ^0 in (|2.1[) . one gets 

det(d t S+V(x,VS)) = 0, 
which leads to the eikonal equation 

d t S ± + h ± (xyS ± )=0. (2.3) 
Applying the projection Il ± (a;,VS' ± ) to (12. 2[) . one gets the solvability condition 

n ± {x,VS ± )(d t + a-V)u^=0. 
After a series of calculations [30l[3T], the following transport equation can be derived: 

d t u± + (w^as.VS*) • V) v£ + \ (V • w ± (a;, V5±)) u± = A ± (x,VS ± )u±, (2.4) 

with 

w ± oc ) *)=v^ ± (s,e), 

k I i 1 

+ ^(*i-A)-((ti-A)-VA + \Vh ± ). 

Remark 2.1. If the external magnetic potential is zero, i.e. A — 0, then 
A ± (x,Z) = -±<x-VV(x) + ^-VV(x). 



3. The Gaussian beam method In this section, we derive the Gaussian 
beam method using both Lagrangian and Eulerian formulations. We first introduce 
the Lagrangian Gaussian beam method for solving the Dirac equation, then discuss 
the Eulerian Gaussian beam method. 

3.1. The Lagrangian formulation In this subsection, we describe how to 
solve the Dirac equation Ijl.ljl by the Lagrangian Gaussian beam method, which is 
given by the following ansatz: 

0f a ± (i ) x,y o ) = M ±(t,y ± )e iT± (^'« ± )/ E ) (3.1) 

where t/ ± =y ± (t,y ) an( i T^^^x^y^) is a second order Taylor truncated phase func- 
tion 



T ± (t,x,y ± ) = S ± (t,y ± )+ti ± (t,y ± )-(x-y ± ) + -(x-y ± ) T M ± (t,y ± )(x-y ± ). 
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Here S ± el, M 3 , M ± € C 3x3 and u± G C 4 . Then the evolutionary ODEs of the 
Lagrangian Gaussian beam (|3.2I) - (I3.5|) can be derived by the eigenvalue decomposition 
technique and the standard Gaussian beam method (for details see Appendix): 

dy ± 



dt 

dt 
dS ± 
dt 
dM ± 
dt 



:V^ ± , (3.2) 

■—Vyh*, (3.3) 

-W € h ± -^ ± -h ± , (3.4) 

: -V yy h ± - V y ^/i ± A/ ± - M ± V^ y /i ± - M ± V^/i* M ± , (3.5) 



^=-l(V y . U) ±) u ±+Au±, (3.6) 

in which ^ y ± =y ± (t,y ),Z ± =Z ± (t,y ± (t,y )),S ± = S ± (t,y ± (t,y )),M ± = 
M ± (t,y ± (t,y )) and Uq =Uq (t,y ± (t,y )). The ODE of Lagrangian Gaussian 
beam amplitude (|3.6j) can be given by the solvability condition and high order 
Gaussian beam formulation. Equations (|3.2l) - (|3.3p are the ray tracing equations, and 
the equation Q3.5P is a Riccati equation which can be alternatively solved by the 
following dynamic first order system 

d 

^ r = (V^ ± )P ± + (V«/ l ± )i? ± , (3.7) 

^ = - (V^) P ± - (V^) R ± . (3.8) 

Then the Hessian matrices satisfy Af ± = i? ± (P ± ) . After that, the Lagrangian 
Gaussian beam solution to the Dirac equation (jl.ll) is constructed as 

*/o(*. aj )= (2^) J R3 ( r e( x - y + ) < >la i^ x ^o) + r e( x - y')4>la (t,x,y )) dy Q , 

(3.9) 

where g C5 (R 3 ) , > is a truncation function with rg = 1 in a ball of radius 9 > 
about the origin. The discrete form of (|3.9[) is given as 



r. 



1 



*f»(*,aO = I ^- J ^r e ( : r-y+(t,y J ))^+(t, a; ,y u )Ay 



+ E r ^-^ + (^o))0L + (*>*>s/u)Ay o , (3.10) 

where y J are the equidistant mesh points, and N yo is the number of the beams 
initially centered at y\ y The initial conditions for the equations (|3.2I) - (I3.6|) are as 
follows QUO!] 

y ± (o,y )=2/o, (3-11) 
| ± (0,y o ) = V^(y o ), (3.12) 

S ± (0,y ) = S , z(y ). (3-13) 
M ± (0,y )=V 2 5 / (y )+ i /, (3.14) 

u±(0,y ) =1^(1/0, V5/( Vo ))tij(yo). (3-15) 



() 
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Remark 3.1. The evolutionary equation p.6l) and the initial condition (|3.15[) ensure 
that 

U + u+(t,y ) = u+(t,y ), U + u (t,y ) = 0, 
n + it(7(ty ) = 0, II~Uo(t,y ) = Uo (t,y ), 

/or Vi>0. The related discussion for the semi-classical limit can be found in \30^ . 

Remark 3.2. As discussed in [10], to compute the Gaussian beam solutions for the 
Schrodinger equation, the optimal mesh can be 0{e^) and the time step requirement 
is of order 0(e^). Therefore, the total computational cost is 0(e~ 2 ) for the 3-D sim- 
ulation. On the other hand, the computational cost for the direct numerical methods 
should be at least 0(e~ 3 ) for accurate physical observables and 0(e~ i ) for accurate 
wave fields. It is obvious to see that the Gaussian beam method is much more efficient 
in the semi- classical regime £<t;l. Since the solutions of the Dirac equation has the 
similar high frequency structures to the ones of the Schrodinger equation. The similar 
discussions can be made to the Dirac equation. 

3.2. The Eulerian formulation In this subsection, the Eulerian Gaussian 
beam method using the level set method 7, 9 is introduced to solve the Dirac equation 
(|1.1[) . The analogous derivations are given in details by the former work of Jin et al 

Define the Liouville operator as 

then the level set equations corresponding to equations (|3.2p - (13.6[) are given by 

£=^±=0, (3.16) 
£ ± S ± =V ( :h ± -€ ± -h ± , (3.17) 

C ± ut = ~{Vyu ± )ut+Au±. (3.18) 

Here <f ± = tp ± (t,y,^) £C 3 are the level set functions. The zero level set of Re[<^] 
gives the (multi- valued) velocity. We also have the phase 5' ± = >S' ± (i,y,£) €R and the 
amplitude — u^(t,y,^) gC 4 in the phase space. To be compatible with the initial 
data ()3.11|) - ()3.15|) . we use the following initial condition: 

<P ± (o,y,£) = -iy+(Z-VyS i (y)), (3.19) 

S ± (0,v,€)=S / (y), (3.20) 
u£(Q,y,t) =n ± ( 2/0 ,V5 / (y ))w 7 (y ). (3.21) 

From p.!6|) and (|3. 191) . the Hessian matrices are constructed via 

M ± = -V y <p(Vscp)-\ 

As a result of this property, we don't need to solve the level set equations for M^, 
or ^ corresponding to equations p.5j) and p.7[) - (|3.8|) as was done in the Eulerian 
Gaussian beam method [TH [T7] . This can save a lot of computational resources 
especially for such 3D problem. After that, the Eulerian Gaussian beam solution to 
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the Dirac equation (jl.ll) is constructed as 

- (g^) 2 J J r e (x~y) (<t>l+(t,x,y,t)6(M<P + (t,y,li)}) 

+ (t> E -(t,x,y,S)6(Re[<p-(t,yS]))dSdy, (3.22) 

where 

T ± (t,x,y,$)=S ± (t,y,£)+Z-(x-y) + ^(x-y) T M ± (t,y,£)(x-y). 

Remark 3.3. The equation (|3.22[) can be solved by a discretized delta function integral 
method \3b\ |#7| / or a local semi-Lagrangian method UOl |i7| /. 

Remark 3.4. One can use the local level set method \21\ \25\/ to solve the level set 
equation in the vicinity of a lower- dimensional zero level curve of i?e[(^ ± ] to reduce 
the total computational cost for the Eulerian Gaussian beam method. An alternative 
efficient way is to use the semi-Eulerian Gaussian beam method proposed in \12^ . 

4. Numerical examples In this section, we present several numerical exam- 
ples to show the accuracy and efficiency of the Gaussian beam method. We compute 
the solution of the Dirac equation (jl.ljl by the time-splitting spectral scheme [6] . The 
reference solutions VP e (t,x) are computed on a very fine mesh and a very small time 
step. In all the numerical examples, the truncation parameter in (|3 . 10[) is chosen 
large enough so that the cut-off error is almost zero. 

Example 1. We consider the zero external fields, i.e. V(x)=0 and A(x)=0. The 
initial condition for the Dirac equation (|l.l|l - (|1.2|l is 

%(x) = e-^ X , X=(1,0,0,0) T , d=^. 
In this example, *| = = *1 = and 

ied t ^l = ^l 

which can be explicitly solved as 

\&f (t,x) = e id* e i . 

The I 1 , 1 2 and errors between the solutions of the Dirac equation , 4' £ and those 
of the Gaussian beam method & e GB for different e are given in Table 14.11 Here we 
take i = 0.5, the time step and mesh size of the Gaussian beam method satisfy At = 
0(52), Ay = 0(ez). We plot the wave amplitudes and absolute errors for different e 
in Figure 14.11 from which, one can see the Gaussian beam methods is more accurate 
for small e and converges nearly first order with respect of e. On the other hand, the 
absolute error could be large for big e, e.g. the absolute l°° error of the Gaussian beam 
solution could be 0.500 for £ = ^q- It is because we are considering the asymptotic 
numerical method. The approximations may not good when the asymptotic parameter 
e is not small enough. However, the l°° error decays almost linearly in e, and we can 
still conclude that the Gaussian beam method is accurate and efficient in the semi- 
classical regime e <C 1. 
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Fig. 4.1. Example 1, at time t = 0.5,£3 = 0, from top to down, they are the amplitude of 

Gaussian beam solutions (Left) \&q b \ and the absolute error (Right) |>t e — for different 

e=^ I 1 — 

256 ' 512 ' 1024 ' 





£ 


1 

2Bfi 


i 

512 


i 

1024 


i 






-*gbIIi 


7.08 x 10~ 2 


4.15 x 10~ 2 


2.27 x 10~ 2 


1.19x10" 


-2 




-*gbII 2 


1.29 x 1CT 1 


7.79 x lO- 2 


4.31 xl0~ 2 


2.30 x 10" 


-2 




-*Gslloo 


5.00 x 10- 1 


3.16x 10" 1 


1.80 x 10" 1 


9.71 x 10" 


-2 



Table 4.1. The I , I 2 and l°° errors of the solutions at t = 0.50 for Example 1. The convergence 
rate in e are 0.8591 in the I 1 norm, 0.8311 in the I 2 norm and 0.7912 in the l°° norm respectively. 



Example 2. We consider the same zero external fields, the initial condition for the 
Dirac equation ()l.lj) - (|1.2[) is 

Sq(x) = ^(l + cos27rxi)(l + cos27rx 2 ), 

X(x) = (~ (y/(d Xl So) 2 + (d X2 S r + 1 + 1) ,0,0, ~ (d Xl S + d X2 S ) 
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Fig. 4.2. Example 2, at time t = 0.38, 23 = 0, from top to down, they are the amplitude of 

Gaussian beam solutions (Left) \&q b \ and the relative error (Right) - — r GB 1 for different e = 

1 1 1 
512 ' 1024 ' 2048 ' 



In this example, ^2 = ^3 = an d 

iedt^i = *i - ied Xl * 4 - ed X2 * 4 , 
i£9 i *4 = -i£9 xl * 1 +ea x2 *i-* 4 , 

which reduces to a two dimensional problem and can be solved by the time-splitting 
spectral method in only one time step. Due to the compressive initial velocity, the 
caustics will form at about issO.56. The I 1 , l 2 and l°° errors between the solutions 
of the Dirac equation and those of the Gaussian beam method & e GB for different 
e are given in Tables |4~2"1|4.3[ for t = 0.375 and t = 0.56 respectively. We remark that 
in Table 14.31 we compare the relative ^°°-error since the caustics form. We plot the 
wave amplitudes and relative errors for different e in Figures 14.2114.31 for which, we 
can draw the same conclusions as in Example 1. 



Example 3 (Harmonic oscillator) . We consider the zero external magnetic poten- 
tial A(x) =0 and the quadratic external electric potential V(x) = ^\x\ . The initial 
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Fig. 4.3. Example 2, at time t = 0.56, xz = 0, from top to down, they are the amplitude of 

Gaussian beam solutions (Left) \&q b \ and the relative error (Right) - — ^ s n GB for different e = 

l l l 

512 ' 1024 ' 2048 ' 
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1 

512 


i 

1024 


i 

2048 


i 

4nsfi 






1.23 X10" 1 


6.41 x 10~ 2 


3.26 x 10~ 2 


1.71 x 10~ 2 






2.21 xlO" 1 


1.23 X10" 1 


6.53 x 10~ 2 


3.37 x 10~ 2 




-*GfllL 


8.36 x 10- 1 


5.40 x 10- 1 


3.09 x 10- 1 


1.57 x lO" 1 



Table 4.2. The I , I 2 andl°° errors of the solutions at t = 0.38 for Example 2. The convergence 
rate in e are 0.9490 in the I 1 norm, 0.9051 in the I 2 norm and 0.8111 in the l°° norm respectively. 



SI 2 



1024 



-*gbII 2 

H* e -*GBjL 



8.62 xlO~ 2 2.92 x 10~ 2 
2.28 xlO" 1 1.26 xlO" 1 
1.76 xlO" 1 1.06 xlO- 1 



i 

2048 



am 

1.14xl0~ 2 5.11 xl0~ 3 
6.62 xlO" 2 3.49 x 10~ 2 
6.14xl0~ 2 3.23 xl0~ 2 



Table 4.3. The I 1 , 1 2 and l°° errors of the solutions at t = 0.56 for Example 2. The convergence 
rate in e are 1.3682 in the I 1 norm, 0.9030 in the I 2 norm and 0.8177 in the i°° norm (relative errors) 
respectively. 
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t=0 t=1 t=2 




Fig. 4.4. Example 3, the amplitude of Gaussian beam solutions at different time, here e = 
gjo, X3 = 0. Note at time t = l and t = 8 the amplitude are cutted since it is too big near caustics. 



condition for the Dirac equation (|l.ll) - (|1.2p is 

(*!-(). l) 2 +(* 2 +0.1) 2 +x§ 1 

¥§(a;) = e 5P X , X=(1,0,0,0) T , d=— . 

This is a full 3D problem. The time splitting spectral method is very expensive 
because of the large requirement of memory for very small e, which the Gaussian 
beam method is accurate. In Figure 14.41 we depict the wave amplitude at different 
time t. In this example, we choose e = . We can see that the wave packet moves 
in circles due to its interaction with the harmonic external potential. 



5. Conclusion In this work, we developed the Gaussian beam method for 
the Dirac equation. The Eulerian Gaussian beam method provides a simple way to 
compute the Hessian matrices for the phase, as in [TU] . The proposed method is shown 
numerically to be accurate and efficient. The required mesh size and time step should 
be of 0(y/e). Compare to the traditional numerical method, which requires mesh 
size to be O(e), the computational cost for our method is cheap, especially when e 
is very small. A more interesting question is to simulate the graphene by using the 
Gaussian beam method. We are currently investigating this important model and 
hope to report our progress in the near future. 
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Appendix. In this appendix, we give the detailed derivation of the Lagrangian 
Gaussian beam method for the Dirac equation. For convenience, we drop the super- 
script ± in (|3.ip 

c/>l(t,x,y )=u(t,x,y)e iT ( t >*<yye, (A-l) 

with 

T{t,x,y) = S(t,y)+Z{t,y)-(x-y) + -(x-y) T M(t,y)(x-y), (A-2) 
u(t,x,y) = u (t,x,y) + eui(t,x,y), (A-3) 

u {t,x,y) = Uo(t,y) + (x-y)-VyU (t,y) + -(x-y) T (Vy y u (t,y)){x-y),(A-4:) 

ui(t,x,y) = ux(t,y). 

Without loss of generality, we can assume that H(x,V y S)uo = UQ. Taking this into 
(jl.ip and matching the first two leading order of e, one obtains 

dy 

dt 

In order to get a nontrivial solution Uq ^ in (IA-5|) , we need 

det (d t T+^--V y T+V{x,V x T)^=0. 
Since h(x,£) is the eigenvalue of V(x,£), one gets 



d t T+-£--V y T + V{xy x T) )«o = 0, (A-5) 
ifatUo + ^-VyUo + a-V.So )-( 9 t T+ V w T+2?(sb,V.T) ) ux = 0. (A-6) 



ftT+^.V w T + /i(x,V 1B T)=0. (A-7) 



Taking the first and second order derivatives with respect to x in (|A-7|) gives 

dv 

d t (V»T) + -jj • V yx T + V x ft + V c ft • V XX T = 0, (A-8) 

5 t (V xx T) + -f--V hV xx T 

+ V^TV^ft + V xx TV^hV xx T + V^ft ■ V XXX T = 0. (A-9) 

Consider (|A-2[) and evaluating (|A-6|) - (|A-9|) at cc = y yield 

^+^-(V„5-€) + ^,€) = 0» (A-10) 

<9t£ + ^ ■ (V„£ - Af ) + V y ft + V 4 ft ■ M = 0, (A-ll) 
dv 

d t M + -| • Vj,M + V yy ft + V y€ ftM + Af V €y ft + MV^hM = 0, (A-12) 



i(9tUo + a-V s uo)+(ft5+^-(V w 5-e)+2?(x,0 )«i = 0. (A-13) 
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We choose the beam center that satisfies 

dy 



d*= V ^ 
then (|A-10P - (|A-12|) can be written as 

— = -Vyyh- V yi hM - MV^yh- MV^hM. 

To obtain the transport equation for uq, we apply the projection H(x,W y S) to (|A-13j) : 

U(x, V y S) (d t + a ■ V y ) m = 0. 

This is the solvability condition for U\ , and one can finally get the ODE for Lagrangian 
Gaussian beam amplitude after a series of calculations 



du 1 
—— = --(\/y-LJ)Uo+Au . 



Remark 5.1. The amplitude (|A-3[) - (|A-4[) is expanded to higher order in both the 
Taylor expansion and the asymptotic expansion, which is different from the Gaussian 
beam method for the Schrodinger equation The reason is that a higher order 

asymptotic expansion is needed for deriving the transport equation for the amplitude 
when using the solvability condition. The higher order Taylor expansion should be 
used to match the high order asymptotic expansion for the Gaussian beam method. 

Remark 5.2. The equation (j A-4|) can be written in a more general form, e.g. 

uo{t,x,y) = UQ (t,y) + (x-y)-u 01 (t,y) + ^(x-y) T Uo2(t,y)(x-y). 

where u m € C°°(R 4 ,C 4 ), Uq\ G C°°(M 4 ,C 3><4 ) and u 02 eC 0O (R 4 ,C( 3x3 ) x4 ). Since 
there are more freedoms than restrictions, one can easily formulate them as 

to close the system. This is consistent to the equation (| A-4[l . 
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